Introduction to neuroimaging data with Python

In this tutorial we will learn how to load, plot, and manipulate neuroimaging data in Python

Written by Luke Chang

Nibabel

Neuroimaging data is often stored in the format of nifti files .nii which can also be compressed using gzip .nii.gz. These files store both 3D and 4D data and also contain structured metadata in the image header.

There is an very nice tool to access nifti data stored on your file system in python called nibabel. If you don't already have nibabel installed on your computer it is easy via pip. First, tell the jupyter cell that you would like to access the unix system outside of the notebook and then install nibabel using pip. You only need to run this once (unless you would like to update the version).

nibabel objects can be initialized by simply pointing to a nifti file even if it is compressed through gzip. First, we will import the nibabel module as nib (short and sweet so that we don't have to type so much when using the tool). I'm also including a path to where the data file is located so that I don't have to constantly type this. It is easy to change this on your own computer.

We will be loading an anatomical image from subject S01 from the open localizer dataset. See this paper for more information about this dataset.

If we want to get more help on how to work with the nibabel data object we can either consult the documentation or add a ?.

The imaging data is stored in either a 3D or 4D numpy array. Just like numpy, it is easy to get the dimensions of the data using shape

We can also directly access the data and plot a single slice using standard matplotlib functions.

Try slicing different dimensions (x,y,z) yourself to get a feel for how the data is represented in this anatomical image.

We can also access data from the image header. Let's assign the header of an image to a variable and print it to view it's contents.

Some of the important information in the header is information about the orientation of the image in space. This can be represented as the affine matrix.

The affine matrix is a way to transform images between spaces.

Here is a short tutorial on affine transformations from the nibabel documentation.

Here is another nice tutorial from nilearn in 2D space.

We have voxel coordinates (in voxel space). We want to get scanner RAS+ coordinates corresponding to the voxel coordinates. We need a coordinate transform to take us from voxel coordinates to scanner RAS+ coordinates.

In general, we have some voxel space coordinate $(i, j, k)$, and we want to generate the reference space coordinate $(x, y, z)$.

Imagine we had solved this, and we had a coordinate transform function $f$ that accepts a voxel coordinate and returns a coordinate in the reference space:

$(x, y, z) = f(i, j, k)$

$f$ accepts a coordinate in the input space and returns a coordinate in the output space. In our case the input space is voxel space and the output space is scanner RAS+.

In theory $f$ could be a complicated non-linear function, but in practice, we know that the scanner collects data on a regular grid. This means that the relationship between $(i, j, k)$ and $(x, y, z)$ is linear (actually affine), and can be encoded with linear (actually affine) transformations comprising translations, rotations and zooms wikipedia linear transform

Scaling (zooming) in three dimensions can be represented by a diagonal 3 by 3 matrix. Here's how to zoom the first dimension by $p$, the second by $q$ and the third by $r$ units:

$$ \begin{bmatrix} x\\ y\\ z \end{bmatrix} \quad = \quad \begin{bmatrix} p & i\\ q & j\\ r & k \end{bmatrix} \quad = \quad \begin{bmatrix} p & 0 & 0 \\ 0 & q & 0 \\ 0 & 0 & r \end{bmatrix} \quad \begin{bmatrix} i\\ j\\ k \end{bmatrix} $$

A rotation in three dimensions can be represented as a 3 by 3 rotation matrix wikipedia rotation matrix. For example, here is a rotation by $\theta$ radians around the third array axis:

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \quad = \quad \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} $$

This is a rotation by $\phi$ radians around the second array axis:

$$ \begin{bmatrix} x \\ y \\ z \\ \end{bmatrix} \quad = \quad \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} $$

A rotation of $\gamma$ radians around the first array axis:

$$ \begin{bmatrix} x\\ y\\ z \end{bmatrix} \quad = \quad \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\gamma) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \cos(\gamma) \\ \end{bmatrix} \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} $$

Zoom and rotation matrices can be combined by matrix multiplication.

Here's a scaling of $p, q, r$ units followed by a rotation of $\theta$ radians around the third axis followed by a rotation of $\phi$ radians around the second axis:

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \quad = \quad \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \quad \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} p & 0 & 0 \\ 0 & q & 0 \\ 0 & 0 & r \\ \end{bmatrix} \quad \begin{bmatrix} i\\ j\\ k\\ \end{bmatrix} $$

This can also be written:

$$ M \quad = \quad \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \quad \begin{bmatrix} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} p & 0 & 0 \\ 0 & q & 0 \\ 0 & 0 & r \\ \end{bmatrix} $$$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \quad = \quad M \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} $$

This might be obvious because the matrix multiplication is the result of applying each transformation in turn on the coordinates output from the previous transformation. Combining the transformations into a single matrix $M$ works because matrix multiplication is associative -- $ABCD = (ABC)D$.

A translation in three dimensions can be represented as a length 3 vector to be added to the length 3 coordinate. For example, a translation of $a$ units on the first axis, $b$ on the second and $c$ on the third might be written as:

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \quad = \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} \quad + \quad \begin{bmatrix} a \\ b \\ c \end{bmatrix} $$

We can write our function $f$ as a combination of matrix multiplication by some 3 by 3 rotation / zoom matrix $M$ followed by addition of a 3 by 1 translation vector $(a, b, c)$

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \quad = \quad M \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} \quad + \quad \begin{bmatrix} a \\ b \\ c \end{bmatrix} $$

We could record the parameters necessary for $f$ as the 3 by 3 matrix, $M$ and the 3 by 1 vector $(a, b, c)$.

In fact, the 4 by 4 image affine array does include exactly this information. If $m_{i,j}$ is the value in row $i$ column $j$ of matrix $M$, then the image affine matrix $A$ is:

$$ A \quad = \quad \begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3} & a \\ m_{2,1} & m_{2,2} & m_{2,3} & b \\ m_{3,1} & m_{3,2} & m_{3,3} & c \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} $$

Why the extra row of $[0, 0, 0, 1]$? We need this row because we have rephrased the combination of rotations / zooms and translations as a transformation in homogenous coordinates (see wikipedia homogenous coordinates). This is a trick that allows us to put the translation part into the same matrix as the rotations / zooms, so that both translations and rotations / zooms can be applied by matrix multiplication. In order to make this work, we have to add an extra 1 to our input and output coordinate vectors:

$$ \begin{bmatrix} x \\ y \\ z \\ 1 \end{bmatrix} \quad = \quad \begin{bmatrix} m_{1,1} & m_{1,2} & m_{1,3} & a \\ m_{2,1} & m_{2,2} & m_{2,3} & b \\ m_{3,1} & m_{3,2} & m_{3,3} & c \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} i \\ j \\ k \\ 1 \end{bmatrix} $$

This results in the same transformation as applying $M$ and $(a, b, c)$ separately. One advantage of encoding transformations this way is that we can combine two sets of [rotations, zooms, translations] by matrix multiplication of the two corresponding affine matrices.

In practice, although it is common to combine 3D transformations using 4 by 4 affine matrices, we usually apply the transformations by breaking up the affine matrix into its component $M$ matrix and $(a, b, c)$ vector and doing:

$$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \quad = \quad M \quad \begin{bmatrix} i \\ j \\ k \end{bmatrix} \quad + \quad \begin{bmatrix} a \\ b \\ c \end{bmatrix} $$

As long as the last row of the 4 by 4 is $[0, 0, 0, 1]$, applying the transformations in this way is mathematically the same as using the full 4 by 4 form, without the inconvenience of adding the extra 1 to our input and output vectors.

You can think of the image affine as a combination of a series of transformations to go from voxel coordinates to mm coordinates in terms of the magnet isocenter. Here is the EPI affine broken down into a series of transformations, with the results shown on the localizer image:

Applying different affine transformations allows us to rotate, reflect, scale, and shear the image.

For example, let's try to reflect the image so that it is facing the opposite direction.

Now maybe we would like to shift the image 10 units in the x direction.

What if we wanted to make the brain smaller by applying a scaling transformation?

How about applying a rotation of 0.3 radians?

Nilearn

There are many useful tools from the nilearn library to help manipulate and visualize neuroimaging data. See their documentation for an example.

Let's make sure it is installed first.

Now let's load a few different plotting functions from their plotting module

Nilearn plotting functions are very flexible and allow us to easily customize our plots

try to get more information how to use the function with ? and try to add different commands to change the plot

nilearn also has a neat interactive viewer called view_img for examining images directly in the notebook.

The view_img function is particularly useful for overlaying statistical maps over an anatomical image so that we can interactively examine where the results are located.

As an example, let's load a mask of the amygdala and try to find where it is located.

We can also load 4D data such as a series of epi images. Here we will load a short functional run from the same particiapnts.

For some reason nibabel is reading this nifti as a 5D image. Let's quickly remove the extra dimension.

Ok, looks like this fixed the dimensions.

Now, let's plot the average voxel signal intensity across the whole brain for each of the 128 TRs.

Notice the slow linear drift over time, where the global signal intensity gradually decreases.

Brainiak

We highly recommend checking out brainiak another great python toolbox developed at Princeton in collaboration with Intel. They also have a number of tutorials for working with advanced data analyses.

nltools

nltools is a python toolbox that makes it easy to work with neuroimaging data. It uses many functions from nilearn and other python libraries. The bulk of the nltools toolbox is built around the Brain_Data() class. The concept behind the class is to have a similar feel to a pandas dataframe, which means that it should feel intuitive to manipulate the data.

The Brain_Data() class has several attributes that may be helpful to know about. First, it stores imaging data in .data as a vectorized features by observations matrix. Each image is an observation and each voxel is a feature. Space is flattened using nifti_masker from nilearn. This object is also stored as an attribute in .nifti_masker to allow transformations from 2D to 4D matrices. In addition, a brain_mask is stored in .mask. Finally, there are attributes to store either class labels for prediction/classification analyses in .Y and design matrices in .X. These are both expected to be pandas DataFrames.

We will give a quick overview of basic Brain_Data operations, but we encourage you to see our documentation for more details.

Brain_Data Basics

To get a feel for Brain_Data let's load an example anatomical overlay image that comes packaged with the toolbox.

To view the attributes of Brain_Data use the vars() function.

Brain_Data has many methods to help manipulate, plot, and analyze imaging data. We can use the dir() function to get a quick list of all of the available methods that can be used on this class.

To learn more about how to use these tools either use the ? function, or look up the function in the api documentation.

Ok, now let's download a more interesting dataset from neurovault.

Here we fetch the pain dataset used in Chang et al., 2015 from neurovault. In this dataset there are 28 subjects with 3 separate beta images reflecting varying intensities of thermal pain (i.e., high, medium, low). The data will be downloaded to ~/nilearn_data, and automatically loaded as a Brain_Data() instance. The image metadata will be stored in data.X.

Here are a few quick basic data operations.

Find number of images in Brain_Data() instance

Find the dimensions of the data (images x voxels)

We can use any type of indexing to slice the data such as integers, lists of integers, slices, or boolean vectors.

Calculate the mean for every voxel over images

Calculate the standard deviation for every voxel over images

Methods can be chained. Here we get the shape of the mean.

Brain_Data instances can be added and subtracted

Brain_Data instances can be manipulated with basic arithmetic operations.

Here we add 10 to every voxel and scale by 2

Brain_Data instances can be copied

Brain_Data instances can be easily converted to nibabel instances, which store the data in a 3D/4D matrix. This is useful for interfacing with other python toolboxes such as nilearn

Brain_Data instances can be concatenated using the append method

Lists of Brain_Data instances can also be concatenated by recasting as a Brain_Data object.

Any Brain_Data object can be written out to a nifti file.

Images within a Brain_Data() instance are iterable. Here we use a list comprehension to calculate the overall mean across all voxels within an image.

Plotting

There are multiple ways to plot data.

Brain_Data() instances can be converted to a nibabel instance and plotted using any plot method such as nilearn.

In addition, there is a fast axial montage .plot() function.

Finally, there is an interactive .iplot() method.

Ok, that's the basics. Brain_Data can do much more!

Check out some of our tutorials for more detailed examples.

We'll be using this tool throughout the course.